VMIproc Runs 128-137

27/11/20 PH

Data analysis for 450/440eV pump-probe, 1\mum focus (I think) - should be comparable to previous analysis posted by TMO team (Taran...?), running from v4 preprocessed data.

Using VMIproc class for VMI data analysis & interactive plotting, including uncertainties.

Implemented:

  • Plotting
  • Inversion
  • Uncertainties (bootstrap protocol)

To do:

  • Masking
  • Covariance

Background:

In [1]:
# # Memory profiler - OPTIONAL for testing
# # https://timothymonteath.com/articles/monitoring_memory_usage/
# %load_ext memory_profiler
# %memit


# # Quick hack for slow internet connection!
# %autosave 36000
In [2]:
# Standard imports
from pathlib import Path
import sys

import numpy as np
import xarray as xr

# Import class from file
sys.path.append(r'/cds/home/p/phockett/dev/tmo-dev/tmo')
import tmoDataBase as tb
# from vmi import VMI  # VMI class
from inversion import VMIproc  # VMI processing class

tb.setPlotDefaults()  # Set plot defaults (optional)

Load data

Set dataset

In [3]:
# Set main data path
baseDir = Path('/reg/data/ana15/tmo/tmolw0618/scratch/preproc/v4')

# Setup class object and which runs to read.
data = VMIproc(fileBase = baseDir, runList = range(128,137+1), fileSchema='_v4')

# The class has a bunch of dictionaries holding info on the data
# data.runs['files']

# Read in the files with h5py
# There is very basic checking implemented currently, just to see if 'energies' is present.
data.readFiles()
*** WARNING: key 131 missing energies data, will be skipped.
Read 10 files.
Good datasets: [128, 129, 130, 132, 133, 134, 135, 136, 137]
Invalid datasets: [131]

Filtering

In [4]:
# Add a filter to an existing filterset
# Set energy filter to match previous analysis
# TODO: implement event count filter here to check for saturation
# data.setFilter({'signal':{'energies':[0.05, 0.1]}})
data.setFilter(reset=True)

data.setFilter({'signal':{'evrs':[1,1,70]}})  # For v4 data, gas == evrs[:,70] true/false
data.setFilter({'bg':{'evrs':[0,0,70]}})
data.setFilter({'energies':[0.05, 0.1]})
data.filters
Out[4]:
{'signal': {'evrs': [1, 1, 70], 'energies': [0.05, 0.1]},
 'bg': {'evrs': [0, 0, 70], 'energies': [0.05, 0.1]}}

Generate VMI images

Note this defaults to electron coords ('xc','yc').

In [5]:
data.genVMIXmulti()

# Check an image with Xarray plotter (Matplotlib)
data.imgStack['signal'].sel(run=128).pipe(np.log10).plot.imshow()
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/cds/home/p/phockett/.conda/envs/ps-4.0.8-local-clone/lib/python3.7/site-packages/xarray/core/computation.py:700: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)
Out[5]:
<matplotlib.image.AxesImage at 0x7f7ec152fc90>
In [30]:
# Check some metrics
print(data.data[129]['signal']['metrics'])
print(data.data[129]['bg']['metrics'])
{'shots': (99398, 1000), 'selected': 85201, 'events': 85201000, 'norm': 99398}
{'shots': (99398, 1000), 'selected': 14197, 'events': 14197000, 'norm': 99398}
In [6]:
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack
Out[6]:
<xarray.Dataset>
Dimensions:  (run: 9, xc: 1048, yc: 1048)
Coordinates:
  * xc       (xc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * yc       (yc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * run      (run) int64 128 129 130 132 133 134 135 136 137
    norm     (run) int64 8733 99398 108119 109603 ... 106583 111458 29521 81092
Data variables:
    signal   (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bg       (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    sub      (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
In [7]:
# View full image stack
data.showImgSet()
WARNING:param.RasterPlot23144: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

Set centre point

This is set as native pixel value, or assumed to be the max point in the image.

In [ ]:
# data.setCentre([553, 546])  # Set explictly
data.setCentre()  # Default case - set as max ind.
data.checkCentre()
In [ ]:
# All the coord parameters are set in `centreInds`
# Note here that `cList` gives the array indicies for the centre point, based on the pixel coords set.
data.centreInds

Basic processing

See the previous demo for info, here we'll just generate subtracted images.

In [ ]:
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack

Image inversion

Multiple methods to be supported... so far only Elio's cpBasex is implemented, and this needs a bit of setup.

In [ ]:
# For cpbasex, set the basis path
basisPath = r'/reg/d/psdm/tmo/tmolw0618/scratch/results/tmo_ana/calc/G_r512_k128_l4.h5'

# Set also local version - pip installed version not working with loadG in testing (version mismatch issue, or something in import routine?)
pbasexPath = r'/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex' #/pbasex'

# Import module + basis functions
data.setMethod(pbasexPath = pbasexPath, basisPath = basisPath)

Invert a set of images by setting filterSet (if not set, all images in self.imgStack will be processed). Optionally, set normalisation parameters for the images, the default is norm={'type':'max', 'scope':'global'}, which normalises to the global maximum (i.e. over all images).

Currently the options are:

type:

  • max
  • sum
  • raw

scope:

  • global
  • local
In [ ]:
# Invert a dataset
data.inv(filterSet='signal')
In [ ]:
# Full results are stored in a dictionay, self.proc[filterSet]
filterSet = 'signal'
data.proc[filterSet].keys()
In [ ]:
# ... where 'xr' contains the results in an Xarray
data.proc[filterSet]['xr']
In [27]:
# For rapid plotting per run, use self.plotSpectra()
data.plotSpectra()
In [28]:
# This can also be set to return a HoloMap for more flexibility
hmap = data.plotSpectra(returnMap = True)
hmap.overlay('BLM').layout('run').cols(1)
Out[28]:

Inversion details

Various other stages and results are also logged to the data strucure...

Images are set to 'filterSet-TYPE', where TYPE is:

  • 'fold': Quadrant filter result if quadrant filter is set. These are useful for checking centre point is OK. (TODO: implement quadrant.unfoldQuadrant for full symmetrized images.)
  • 'fit': Reconstructed/fitted images.
  • 'inv': Inverted images

These can be accessed and plotted with the usual functionality.

NOTE: there is currently no default radial masking set, so the default plot scaling will likely be off - adjust with the histogram.

In [ ]:
# Quad filter 
data.showImg(name='signal-fold')
In [ ]:
data.showImg(name='signal-inv')
In [ ]:
# Raw cpbasex output data
data.proc['signal']['pbasex'].keys()

Uncertainties

A basis Poissonian bootstrapping routine is implemented for estimating uncertainties. This operates on the event data, and allows error propagation through the image generation & analysis, but is (consequently) a little slow. (For more on the method employed, see the wikipedeia page#Poisson_bootstrap).)

Note this routine is a little crude at the moment! It will run image generation for all filterSets, run image subtraction, then image inversion & stats for all filterSets, or just the listed element.

In [20]:
# Run bootstrapping & inversion routine, 5 iterations, and calculate stats for the subtracted data only.
data.bootstrapInv(N = 5, filterSet = 'sub')
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 2 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 3 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 4 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 5 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
In [21]:
# Stats are output to a dataset, with variables per iteration
data.stats
Out[21]:
<xarray.Dataset>
Dimensions:  (BLM: 3, E: 512, run: 9)
Coordinates:
  * BLM      (BLM) int64 0 2 4
  * E        (E) float64 0.0 0.000359 0.001436 0.003231 ... 93.01 93.38 93.74
  * run      (run) int64 128 129 130 132 133 134 135 136 137
    XS       (E, run) float64 0.0 0.0 0.0 0.0 0.0 ... 3.243 5.57 1.329 4.304
    mask     (E, run) bool False False False False False ... True True True True
    norm     (run) float64 124.3 124.4 190.2 211.6 115.4 117.8 186.2 43.7 157.7
Data variables:
    0        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6191 -0.5535 -0.5792
    1        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.563 -0.4151 -0.5637
    2        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.5544 -0.4514 -0.6345
    3        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.5808 -0.5251 -0.6106
    4        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.5982 -0.4166 -0.5915
In [22]:
# Similarly to plotSpectra(), there is a plotUncertainties() function for this data
data.plotUncertainties()
In [23]:
# As before, an object can be returned for further plotting options
hmap = data.plotUncertainties(returnMap=True)
hmap.overlay('BLM') #.cols(1)
Out[23]:
In [24]:
hmap.overlay('run').layout('BLM').cols(1)
Out[24]:

Versions

In [25]:
# tmo-dev class version
data.__version__
Out[25]:
'0.0.1'
In [26]:
import scooby
scooby.Report(additional=['holoviews', 'xarray'])
Out[26]:
Fri Nov 27 13:39:51 2020 PST
OS Linux CPU(s) 16 Machine x86_64
Architecture 64bit RAM 125.6 GB Environment Jupyter
Python 3.7.8 | packaged by conda-forge | (default, Jul 31 2020, 02:25:08) [GCC 7.5.0]
holoviews 1.13.5 xarray 0.16.1 numpy 1.19.4
scipy 1.5.3 IPython 7.19.0 matplotlib 3.3.2
scooby 0.5.6